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Abstract 

Background: Haldane's Rule, the tendency for the heterogametic sex to show reduced fertility in hybrid crosses, can 
obscure the signal of gene flow in mtDNA between species where females are heterogametic. Therefore, it is important 
when studying speciation and species limits in female-heterogametic species like birds to assess the signature of gene 
flow in the nuclear genome as well. We studied introgression of microsatellites and mtDNA across a secondary contact 
zone between coastal and interior lineages of Western Scrub-Jays (Aphelocoma califomica) to test for a signature of 
Haldane's Rule: a narrower dine of introgression in mtDNA compared to nuclear markers. 

Results: Our initial phylogeographic analysis revealed that there is only one major area of contact between coastal and 
interior lineages and identified five genetic clusters with strong spatial structuring: Pacific Slope, Interior US, Edwards 
Plateau (Texas), Northern Mexico, and Southern Mexico. Consistent with predictions from Haldane's Rule, mtDNA showed 
a narrower dine than nuclear markers across a transect through the hybrid zone. This result is not being driven by 
female-biased dispersal because neutral diffusion analysis, which included estimates of sex-specific dispersal rates, also 
showed less diffusion of mtDNA. Lineage-specific plumage traits were associated with nuclear genetic profiles for 
individuals in the hybrid zone, indicating that these differences are under genetic control. 

Conclusions: This study adds to a growing list of studies that support predictions of Haldane's Rule using dine analysis of 
multiple loci of differing inheritance modes, although alternate hypotheses like selection on different mtDNA types 
cannot be ruled out. That Haldane's Rule appears to be operating in this system suggests a measure of reproductive 
isolation between the Pacific Slope and interior lineages. Based on a variety of evidence from the phenotype, ecology, 
and genetics, we recommend elevating three lineages to species level: A califomica (Pacific Slope); A. woodhouseii (Interior 
US plus Edwards Plateau plus Northern Mexico); A. sumichrasti (Southern Mexico). The distinctive Edwards Plateau 
population in Texas, which was monophyletic in mtDNA except for one individual, should be studied in greater detail 
given habitat threat. 
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Background 

The study of speciation has been influenced during the 
last decade by increasing awareness that different parts of 
the genome tell different evolutionary stories [1-3]. Previ- 
ously, speciation studies relied heavily on mitochondrial 
DNA (mtDNA), given the numerous benefits offered by 
these data for inferring evolutionary history at diverse 
timescales (e.g., high substitution rates) [4]. Yet, along with 
these benefits come drawbacks [5], among them that 
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mtDNA is a single recombinational unit that can present 
a biased snapshot of evolutionary history when mtDNA 
data are not balanced by independent evidence from nu- 
clear markers. 

One case in which mtDNA data can be particularly 
misleading occurs during gene flow between popula- 
tions in secondary contact. Haldane was the first to note 
that when two differentiated lineages come back into re- 
productive contact, the heterogametic sex (i.e., the one 
with two types of sex chromosomes) often shows reduced 
fertility [6]. The mechanistic reasons behind Haldane's 
Rule, as it came to be known, are under intense study 
[7-10]. While we may not understand why Haldane's Rule 
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occurs, the consequences for female-heterogametic species 
like birds and butterflies are clear: if female hybrids show 
reduced fertility, then mtDNA, which is inherited as a sin- 
gle recombinational unit from the female parent, is less 
likely to be passed on to the next generation. As a result, 
mtDNA might show little or no evidence of hybridization, 
despite the fact that hybridization occurs and paternal, nu- 
clear alleles move between lineages. 

The potential mismatch between gene flow estimated 
from mtDNA and gene flow estimated using nuclear 
markers is not just theoretical; it has been demonstrated 
empirically by studies of female-heterogametic organisms 
[11-18]. This phenomenon is of interest to biologists, 
among other reasons, because of the continuing influence 
of the Biological Species Concept (BSC) and its emphasis 
on cessation of gene flow as a necessary criterion for spe- 
cies delimitation [19]. While the BSC may not be applied 
as commonly in practice as is claimed [20], it is neverthe- 
less evident that many attempts to split species taxonom- 
ically fail due to the existence of gene flow and the 
perceived importance of the BSC. As empirical evidence 
continues to mount demonstrating mito-nuclear discord- 
ance [21,22], the use of nuclear markers has become es- 
sential to species delimitation and the study of speciation. 

Western Scrub-Jays (Aphelocoma californica) are an ex- 
cellent model for testing the mismatch between mtDNA 
and nuclear markers across a hybrid zone predicted by 
Haldane's Rule because previous work documented a zone 
of secondary contact and hybridization between two diver- 
gent lineages. Western Scrub-Jays are widely distributed 
in North America from British Columbia, Canada to 
Oaxaca, Mexico and east to the Great Plains (Figure 1). 
Prior research delimited several lineages based on 
morphology, allozymes, and mtDNA variation [23-26]. 
Two of these lineages - one found primarily in oak 
woodlands along the Pacific slope and the other found in 
pinyon-juniper or pine-oak habitats in the interior US and 
Mexico - meet and hybridize in a few mountain ranges in 
western Nevada east of Lake Tahoe (Figure 2), but are 
otherwise separated from one another by the high- 
elevation conifer forests of the Sierra Nevada. 

These two lineages ('coastal' and 'interior') are thought 
to have diverged from each other anywhere from 1 to 4 
million years ago (Ma) based on divergence times cali- 
brated with a fossil and estimated molecular substitution 
rates [24] . Coastal individuals are easily distinguished from 
those in the interior by their brighter blue plumage, bolder 
white supercilium, and more pronounced blue collar. 
Coastal individuals also have stout, hooked bills, which are 
better for husking acorns than the smaller, more pointed 
bills of interior individuals, which are better suited for 
extracting pine seeds from cones [27,28] . 

The hybrid zone was originally described by Pitelka 
[25] who traversed the zone and analyzed numerous 



museum specimens. He found intermingling of coastal 
and interior phenotypes in western Nevada, as well as in- 
dividuals that appeared to be Fl hybrids. The hybrid zone 
has never been studied in detail using genetic techniques, 
although previous work documented a relatively rapid 
turnover of mtDNA types and mismatches of taxonomy 
and genotype of birds from this zone [23,24]. Neverthe- 
less, a prior proposal to elevate the coastal and interior lin- 
eages to species status [29] failed because of doubts 
concerning the extent of gene flow in the contact zone. 

The purpose of this study was to test for a prediction of 
Haldane's Rule that there will be a mismatch between 
mtDNA and nuclear gene flow across the contact zone. 
Specifically, we predict that if Haldane's Rule is operating 
in this system, we will see a narrower (i.e., steeper) cline in 
the transition of mtDNA haplotypes compared to nuclear 
markers. Before addressing this hypothesis, however, we 
needed to determine the exact nature of the contact zone 
in western Nevada and whether there were any unde- 
scribed contact zones nearby that might confound our re- 
sults. We therefore first assessed broad-scale patterns of 
divergence and gene flow among Western Scrub-Jays 
across their entire geographic range, using deep sampling 
of populations and a suite of nuclear markers. Previous 
molecular studies sampled only one or a few individuals of 
different populations, focusing largely on allozymes [26] 
or mtDNA [23]. Our study design includes nearly 700 in- 
dividuals sampled from throughout the range of the spe- 
cies, thus providing a comprehensive portrait of genetic 
divergence and allowing us to locate potentially unde- 
scribed contact zones. Finally, we correlated proportion 
of hybrid ancestry with phenotypic traits for a subset of in- 
dividuals in the contact zone, as a way to control for the 
influence of environment and therefore test for a genetic 
basis of phenotypic differences between the coastal and in- 
terior lineages. If phenotypic traits are genetically based, 
we expect to see a correlation between neutral genetic as- 
signment to coastal or interior lineages and the phenotypic 
differences known to characterize those lineages. 

Results 

Geographic structure of coastal and interior mtDNA types 

To assess gene flow throughout the range of Western 
Scrub-Jays, we conducted simple mtDNA typing of all 689 
individuals using restriction digest, which showed that 
most populations are composed entirely of coastal or in- 
terior mtDNA types, with very few populations showing 
mixing of mtDNA types. The four populations showing 
mixing were (1) the previously described contact zone in 
western Nevada, which featured several populations hav- 
ing mixed types 100 km east and southeast of Reno and 
Carson City (sites inside the inset square in Figure 1); (2) a 
single interior type (among n = 16 coastal types) found at 
the base of the western Sierra 100 km east of Stockton, 
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Figure 1 Sampling locations for Western Scrub-Jays. Green shading represents heavily forested regions that correspond with high-elevation 
barriers to dispersal. Numbers correspond to site numbers listed in Table 2. Dot color corresponds to the five nuclear genetic clusters from the 
Structure results shown in Figure 3. The square outline denotes the region of the hybrid zone depicted in greater detail in Figure 2. Dots with 
black points inside them represent the three coastal sites away from the hybrid zone that contain one or more interior mtDNA haplotypes. The 
scrub-jay on the lower left is representative of the coastal type, with bolder blue plumage, a larger blue collar, a bolder eyestripe, whiter undertai 
coverts, and a larger, more hooked bill compared to the interior type shown at upper right. Note that the range of the species extends into 
southwestern Canada and into some parts of the Great Plains that are not shown. 



CA (site 6); (3) an even mix of coastal and interior types 
in the eastern Sierra Nevada to the west of the Owens 
Valley near Independence, CA (site 18); and (4) one re- 
markable interior type (among n = 16 coastal types) in 
the Santa Ana Mountains in Orange County, Los 
Angeles (site 23). 

Broad-scale geographic structure of microsatellite 
variation 

To complement the mtDNA typing results with infor- 
mation from nuclear markers, we genotyped all 689 in- 
dividuals at 14 microsatellite loci. No pairs of loci were 
physically linked based on linkage tests. Of the 14 loci, 
one locus (MJG8) showed significant deviations from 
HWE across many populations following Bonferroni 
correction. We removed this locus from the data set. 
Four other loci showed deviations from HWE across 
some populations. We ran downstream analyses with 



and without these loci to investigate whether they were 
driving any of the patterns of genetic structure we observed. 
Results among analyses including and excluding loci were 
qualitatively similar, so we retained the larger group of 13 
loci. No loci showed consistent signs of other microsatellite 
artifacts like large-allele drop-out or null alleles. 

Initial Structure runs from K = 1-30 found increasing 
likelihood (LnL) values until K = 9, followed by a plateau in 
LnL until K = 13, followed by a decrease from K = 14-30. 
Structure runs on successively smaller clusters revealed 
five distinct geographic clusters of nuclear DNA variation 
(Figure 3) that were largely uniform in their population as- 
signment. These genetic clusters were (1) a Pacific Slope 
group; (2) an Interior US group; (3) a group from the 
Edwards Plateau in Texas; (4) an interior group from 
Mexico north of the Transvolcanic Belt (Northern Mexico); 
and (5) an interior group from Mexico south of the Trans- 
volcanic Belt (Southern Mexico). 
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Figure 2 Sampling sites in the hybrid zone. Green shading 
corresponds to heavily forested regions considered to be a barrier to 
dispersal. Numbers correspond to site numbers listed in Table 2. The 
shading of the dots represent their cumulative Structure population 
assignment (including only individuals along the transect run at K=2) 
averaged across all individuals in each population. 



Phylogeny from mtDNA 

To determine if the novel genetic clusters found with 
microsatellites (e.g., the Edwards Plateau group) could also 
be found in a mtDNA phylogeny, we sequenced a subset of 
individuals at the cytochrome b gene (cyt b), including sam- 
pling from throughout the genus. We chose samples to 
both capture the major genetic lineages within Western 
Scrub-Jays and outgroups, and also to sample deeply within 
certain microsatellite groups of interest (e.g., we sequenced 
all individuals from the Edwards Plateau). The larger phyl- 
ogeny, including all Aphelocoma samples, is available 
through Dryad [30]. In Figure 4, we show a trimmed ver- 
sion of this phylogeny including only scrub-jays, which con- 
firms that the current concept of the Western Scrub-Jay 
is paraphyletic because it includes the Island Scrub-Jay 
(A insularis). In addition to the Island Scrub-Jay, phylogen- 
etic analyses placed three other highly supported and diver- 
gent clades within Western Scrub-Jays: (1) a Pacific Slope 
clade that was sister to the Island Scrub-Jay (1.0 posterior 
probability or PP); (2) a clade from Mexico south of the 
Transvolcanic Belt (Southern Mexico; 1.0 PP); (3) a clade 
of all other interior individuals (0.94 PP). Within this last 
clade, we recovered a weakly supported (0.87 PP) group 
from the Edwards Plateau in Texas, with the exception of 
one individual that grouped elsewhere in the interior clade 
(small arrow by terminal tip in Figure 4). There was no 
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Figure 3 Results from Structure runs on 13 nuclear microsatellite loci for successively smaller genetic clusters (K=2 for each run). 

Numbers correspond to sites listed in Table 2. (A) The first run split the individuals broadly into coastal and interior groups, with mixed 
assignment for individuals in the contact zone, which were removed from further runs. (B) The coastal group showed little geographic structure. 
(C) The interior group split into a US Interior cluster and a Mexico cluster. (D) An Edwards Plateau group split from the rest of US Interior; (E) A 
northern Mexico group split from a southern Mexico group; (F) The southern Mexico group showed some evidence for differential assignment 
between Oaxaca and Guerrero populations. 
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Figure 4 BEAST consensus tree of scrub-jays based on cyt b mtDNA sequences new to this study combined with existing sequences 
from GenBank. A more comprehensive version of this tree including all outgroup samples is available through Dryad [30]. The Western Scrub-Jay is 
paraphyletic because it does not include the Island Scrub-Jay. The one texana individual that does not group with the clade of texana individuals with 
0.87 PP is denoted with a small arrow next to the tip label. 



support in mtDNA for a split between Interior US and 
Northern Mexico (as was found in microsatellites), al- 
though some geographic localities grouped together 
with high PP. 



Geographic cline analysis 

Both the mtDNA and microsatellite data showed a step 
cline pattern of variation across the hybrid zone, and the 
null model of no clinal variation had much higher AICc 
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scores in each case (mtDNA: null model = 96.11; cline 
model = 13.76; nuclear DNA: null model = 63.53; cline 
model = 15.72). Model 1 with no scaling was accepted as 
the most likely model for both datasets. Estimates of the 
center of the cline (c) were similar and had overlapping 
confidence intervals. The center of the mtDNA cline 
was ~300 km (280 - 319), and the center of the nuclear 
DNA cline was -315 km (271 - 377). The width of the 
mtDNA cline was -131 km (76 - 270), and the width of 
the nuclear DNA cline was -331 km (145 - 678). The 
estimated width of the mtDNA cline fell outside the 2 
log likelihood confidence intervals for the nuclear DNA, 
and vice versa, therefore the mtDNA cline was significantly 
narrower than the microsatellite cline (Figure 5). Disparity 
in dispersal distances between males and females did not 
account for the difference in cline width between mtDNA 
and microsatellites because our time-of-contact estimates 
under a model of neutral diffusion, which included esti- 
mates of male and female dispersal distance (females: o~ = 
1.5 - 5.0 km; males: o = 0.5 - 3.0 km; [31]), did not overlap 
(mtDNA = -690 - 7600 years, microsatellites = -12,000 - 
440,000 years). In other words, even taking sex-biased dis- 
persal into account, mtDNA showed less diffusion. 

Correlation between phenotype and genotype 

We observed strong correlations between hybrid values 
(Q-scores) from Structure and many phenotypic traits for 
the full data set of 123 individuals that included many pure 
coastal and interior individuals outside the hybrid zone, es- 
pecially for plumage traits and principal components (PC) 
axes describing bill shape (Table 1; see Additional file 1 for 
PCA loadings table). However, by including individuals 
from allopatric ranges, this test potentially confounds gen- 
etic and environmental causes of trait variation. Within the 
hybrid transect and inside the core hybrid zone, where indi- 
viduals experience the same environment, plumage traits, 
but not morphological traits, were correlated with Q-scores 
(Table 1). 

Discussion 

Haldane's Rule and the mismatch in cline width between 
mtDNA and microsatellites 

Our results support a key prediction from Haldane's 
Rule that mtDNA will show a steeper cline across a hy- 
brid zone than nuclear markers. Although Haldane's 
Rule may explain the clinal pattern we observed, this 
pattern could also result from divergent selection on 
mtDNA types [32] or sex-biased dispersal, wherein fe- 
males (and the mitochondrial genomes they pass on) 
tend to disperse less than males. Although we have not 
looked at selection, what little is known about natal dis- 
persal in Western Scrub-Jays suggests that, to the con- 
trary, females disperse about three times farther than 
males on average [31]. The results from our time-since- 
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Figure 5 Transition in genetic types across the hybrid zone for 
(A) mtDNA and (B) nuclear markers. Black circles represent 
sampling localities, the dark line is the cline estimate and the gray 
shading represents confidence around the cline estimate. 



contact analyses, which incorporated estimates of dispersal 
distances, also support predictions from Haldane's Rule, with 
less diffusion of mtDNA compared to nuclear markers. 

Our results showing steeper clines in mtDNA compared 
to nuclear markers add to a gathering corpus of research 
showing discordant rates of DNA transmission across hy- 
brid zones among different marker types [11-18,22], a pat- 
tern that has often been attributed to Haldane's Rule [21]. 
Whatever the root cause, these results caution that con- 
clusions about individual-level gene flow made on the 
basis of mtDNA alone may be underestimated. In our 
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Table 1 Regression (R 2 ) values between hybrid index and 



phenotypic traits within and outside the hybrid zone 


Trait 


All 


Transect 


Core hybrid 


Wing 


0.04* 


ns 


ns 


Tail 


ns 


ns 


ns 


Tarsus 


ns 


ns 


ns 


Bill length 


ns 


ns 


ns 


Bill depth 


ns 


ns 


ns 


Bill width 


0.11*** 


ns 


ns 


Morph PCI ' 


ns 


ns 


ns 


Morph PC2 2 


0.06** 


ns 


ns 


Morph PC3 3 


0.05** 


ns 


ns 


Collar 


0.38*** 


0.23** 


ns 


Vent 


0.51*** 


0.46*** 


0.36*** 


Eyestripe 


0.31*** 


0.25*** 


0.37*** 



1 PC1 describes general size. 

2 PC2 describes a trade-off between body and bill size. 

3 PC3 describes a trade-off between bill length and bill width. 

ns = not significant, *p<0.05, **p<0.01, ***p < 0.001. 



study, gene flow is generally restricted in both mtDNA 
and nuclear markers, due to the low dispersal distances of 
Western Scrub-Jays. In other species, such as long- 
distance dispersers with no sex-biased dispersal, mis- 
matches could be much more pronounced. 

Recent work on the genetic basis for Haldane's Rule sug- 
gests it results from incompatibilities among genes that 
are exposed to selection on hybrid female Z chromo- 
somes, causing sterility [7], whereas only those incompati- 
bilities that are dominant are visible to selection in hybrid 
ZZ males. If Haldane's Rule is indeed the cause of the cli- 
nal mismatch we observed, then our results suggest that 
the coastal and interior lineages of Western Scrub-Jay 
have begun the process of reproductive isolation. Ultim- 
ately, Haldane's Rule as the cause of the mismatch in 
marker introgression rates could be determined through 
field studies or captive mating experiments that document 
the reproductive viability of female offspring resulting 
from hybrid crosses [33]. 

Clarifying the hybrid zone between coastal and interior 
Western Scrub-Jays 

Our results also help clarify the nature of the hybrid zone 
between coastal and interior lineages of Western Scrub- 
Jays. The fact that our hybrid transect populations could 
be fit to a cline suggests that Pitelka [25] was correct that 
the major axis of gene flow is from the northwest to 
southeast, as coastal populations continue around the 
north side of the Sierra Nevada and meet with interior 
populations in western Nevada. However, the coastal in- 
fluence of the Woodford site (45 in Figure 2) warrants fur- 
ther study. In addition to gene flow from the north, some 
gene flow may link the southern end of the Sierra Nevada 



to the Woodford site, with coastal individuals inhabiting 
canyons of the southeastern Sierra that were not sampled 
for this study. 

Our results also support Pitelka's [25] conclusions that 
the hybrid zone is narrow and the result of secondary con- 
tact. Even the oldest time estimates from our time-since- 
contact analysis across the hybrid zone (~437,000 years) 
are younger than divergence estimates of the original split 
between coastal and interior lineages [1 to 4 Ma; 24], sup- 
porting the hypothesis that the hybrid zone results from 
secondary contact and not in situ divergence. In terms of 
broader biogeographic evidence, several other species have 
secondary contact zones where coastal and interior line- 
ages meet in this general region of the Great Basin [34,35]. 

Our results also confirm that the contact zone in west- 
ern Nevada is the only place where coastal and interior 
lineages come into contact to any great degree. Outside of 
western Nevada, two sites (6 and 23 in Figure 1) appear to 
contain examples of rare, long-distance dispersers from 
pure interior populations because the nuclear DNA profile 
of the two individuals matches their mtDNA type (the in- 
terior mtDNA individual at site 6 was 97% interior in nu- 
clear DNA, and the interior mtDNA individual at site 23 
was 95% interior in nuclear DNA). In contrast, at a third 
site in the eastern Sierra Nevada (18 in Figure 1), the three 
individuals with interior mtDNA had strongly coastal nu- 
clear DNA profiles (91%, 95%, and 89% coastal), suggest- 
ing that they were advanced backcrosses from older 
hybridization events (i.e., nuclear DNA profiles of Fls 
would be closer to 50/50). Thus, although there appears to 
be one major contact zone, the exceptions show that the 
reality of dispersal and gene flow among jays in and 
around the Sierras - and across the Mojave Desert farther 
south - is more complex. Fine-scale sampling, combined 
with habitat-based methods of identifying dispersal corri- 
dors [36], could provide a more complete picture of gene 
flow across and around the Sierra Nevada for this and 
other species. 

Previous studies documented phenotypic differences be- 
tween coastal and interior Western Scrub-Jays [25,28], but 
ours is the first to provide evidence for a genetic basis for 
some traits. There was a correlation between nuclear gen- 
etic variation and phenotypic traits, such as bill shape and 
plumage traits, when considering all individuals, even 
those from pure populations (Table 1). However, including 
individuals from different geographic locations potentially 
confounds environmental and genetic sources of pheno- 
typic variation. Still, many of these relationships remained, 
for the plumage traits at least, within the contact zone 
where there are hybrid individuals and all variety of back- 
crosses. The correlation between proportion of interior vs. 
coastal genetic ancestry and phenotypic traits for individ- 
uals that were fledged in the hybrid zone and therefore ex- 
perienced the same environment since fledging provides 
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evidence for a genetic basis of plumage traits, at the very 
least. Why we did not observe similar relationships for bill 
traits is not clear, but could be due to low statistical power, 
stabilizing selection, or greater plasticity in bill traits. Fur- 
ther studies with more individuals and better genomic 
coverage will help discriminate between these alternate 
hypotheses, as will studies that attempt to identify the gen- 
etic loci or environmental pressures responsible for the 
phenotypic differences. 

Genetic differentiation across North America's 
scrub woodland 

Previous genetic studies of species inhabiting North 
America's pine-oak and scrub woodland have uncovered 
unexpectedly deep phylogenetic structure within species 
[37,38]. Results of the current study, which are the most 
comprehensive in terms of geographic, individual, and 
genetic sampling of Western Scrub-Jays, indicate at least 
five genetic clusters with strong spatial structuring (col- 
ored dots in Figure 1). The Pacific Coast and Southern 
Mexico groups have been recovered previously in stud- 
ies of phenotype, allozymes, and mtDNA [23-26]. The 
Edwards Plateau group and the split between Interior 
US and Northern Mexico individuals are novel, although 
the inferred population clusters align with current subspe- 
cies designations [25]. 

This study is the first to show clear evidence of nu- 
clear and mtDNA differentiation of Western Scrub-Jays 
on the Edwards Plateau in Texas (currently a separate 
subspecies, A. c. texana). All individuals in the Edwards 
Plateau population except one possess a unique mtDNA 
haplotype differing from other interior US haplotypes by 
two substitutions. Nuclear and mtDNA differentiation on 
the Edwards Plateau is interesting given the history of isola- 
tion of this region and its unique flora and fauna. The oak 
savanna and Ashe juniper woodland of the Edwards Plateau 
are separated from other scrub woodland to the west by an 
extension of the Chihuahuan Desert in the Pecos Valley. 
Fossil and pollen data suggest, however, that pygmy wood- 
land was continuous across this region, and connected the 
Edwards Plateau to western forests throughout the last gla- 
cial maximum until about 11,000 years ago [39]. At least 
one other distinct bird species with a similar geographic 
range occurs on the Edwards Plateau, the Golden-cheeked 
Warbler (Setophaga chrysoparia), and other organisms in 
this region have distinct populations or subspecies [40,41]. 
Our results lend additional support to the distinctiveness of 
the biota in this region, where native habitat is under per- 
sistent conservation threat from urbanization. 

The nuclear genetic break between the Interior US and 
Northern Mexico populations, occurring roughly where 
the Sierra Madre Occidental and Oriental turn into a net- 
work of sky islands near the US border, could be explained 
either by a break in suitable habitat or a sampling artifact. 



Addressing the latter point, observations posted on eBird 
(ebird.org) confirm observations of Western Scrub-Jays in 
southeastern Arizona and northern Sonora (in the vicinity 
of the genetic break), but these records could represent 
wandering individuals. Likewise few observations exist in 
the vicinity of the genetic break in southwestern Texas, and 
virtually no observations exist from northern Coahuila, 
where congeneric Mexican Jays (Aphelocoma wollweberi) 
occur at lower elevations, seemingly occupying the Western 
Scrub-Jay niche [25,42-44]. Thus, while available evidence 
suggests that the genetic break between Interior US and 
Mexico populations relates to a real distributional break 
near the US border, the exact causes of this break remain 
unclear. Future niche modeling could be used to determine 
if there are cryptic breaks in suitable habitat. 

Lack of genetic differentiation along the Pacific Coast 

In contrast to the highly structured genetic portrait of 
the interior US and Mexico, Western Scrub-Jays along 
the Pacific Slope show a remarkable lack of genetic struc- 
ture despite 2,500 km of nearly unbroken sampling from 
northern Oregon to the southern tip of Baja California. 
Lack of genetic structure over large distances suggests two 
possible scenarios: (1) high levels of dispersal and gene flow 
link coastal populations; or (2) recently expanded coastal 
populations have had insufficient time to accumulate gen- 
etic differences. Available evidence suggests that Western 
Scrub-Jays generally have low dispersal [31], notwithstand- 
ing a few rare, long-distance dispersal events discussed 
above. Furthermore, the original split between coastal and 
interior lineages occurred between 1-4 million years ago 
[24], providing ample time for genetic structure to 
emerge in coastal populations, absent other demographic 
events. Thus, the hypothesis that seems most likely is a re- 
cent expansion following a population bottleneck. How- 
ever, Pleistocene niche models of Western Scrub-Jay 
niches suggest ample habitat along the Pacific Slope dur- 
ing the last glacial maximum [45]. Further study is needed 
to determine the cause of the apparent genetic homogen- 
eity of Pacific Slope populations. 

Taxonomic recommendations 

We recommend splitting the current concept of the Western 
Scrub-Jay into three species, as has already been suggested 
previously [46]. The Pacific Slope lineage (A calif ornica, in- 
cluding subspecies californica, oocleptica, caurina, obscura, 
hypoleuca, superciliosa, immanis, and cactophila) is pheno- 
typically distinct [25], possesses a bill morphology adapted to 
local resources [27,28], is monophyletic in mtDNA, and is 
sister to another recognized species, the Island Scrub-Jay 
[23,24]. This study shows that the Pacific Slope lineage is also 
well-differentiated in nuclear markers and has only a narrow 
zone of contact with interior populations. The genetic clines 
across the hybrid zone are generally steep, and the mismatch 
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between mtDNA and nuclear clines suggests a measure of 
reproductive isolation. These results would appear to meet 
criteria for species distinction under even the most stringent 
species concept. 

Southern Mexico populations (A. sumichrasti, includ- 
ing subspecies sumichrasti and remota) also meet criteria 
for species recognition. They are phenotypically distinct 
[25], with behavioral evidence that some southern pop- 
ulations are cooperative breeders [47], whereas other 
Western Scrub-Jays breed in territorial pairs. This study 
demonstrates that these southern populations are recip- 
rocally monophyletic to other interior populations in 
mtDNA and well-differentiated in nuclear markers. 

Finally, we recommend that the remaining interior popu- 
lations (A. woodhouseii including subspecies texana, wood- 
houseii, nevadae, grisea, and cyanotis) be elevated to species 
level, given their diagnostic phenotype [25], adaptive bill 
morphology [28], high support for monophyly in mtDNA in 
previous work [23,24] and this study (0.94 PP), and differen- 
tiation in nuclear markers (this study). Within this group, 
we do not recommend elevating the genetically distinctive 
Edwards Plateau population to species level at this time be- 
cause our study sampled only one population. More infor- 
mation is needed on potential gene flow with other interior 
US populations, as well as a modern multivariate analysis of 
their phenotypic differentiation and studies of their eco- 
logical and behavioral divergence. Likewise, we do not cur- 
rently recommend splitting Interior US from Northern 
Mexico populations, unless, at minimum, further sampling 
near the US-Mexico confirms that the break seen in micro- 
satellite variation is not the result of a sampling artifact. 

Conclusions 

We report a mismatch in cline widths between mtDNA 
and nuclear markers across a hybrid zone in Western 
Scrub-Jays. This result is consistent with predictions based 
on Haldane's Rule, where there is less mtDNA introgression 
between lineages in female-heterogametic species. This 
study adds to a growing list of studies that support predic- 
tions of Haldane's Rule using cline analysis of multiple loci 
of differing inheritance modes. Structure results of micro- 
satellite variation confirm that there is only one major area 
of contact between coastal and interior lineages. Further- 
more, Structure identified at least five clusters of individuals 
with strong spatial structuring: Pacific Slope, Interior US, 
Edwards Plateau, Northern Mexico, and Southern Mexico. 
Based on a variety of evidence from the phenotype, ecology, 
and genetics, we recommend elevating three lineages to 
species level: A. californica (Pacific Slope); A. woodhouseii 
(Interior US plus Edwards Plateau plus Northern Mexico); 
A. sumichrasti (Southern Mexico). The distinctive Edwards 
Plateau population in Texas, which was also monophyletic 
in mtDNA except for one individual, should be studied in 
greater detail given conservation concerns. 



Methods 

Specimens and DNA extraction 

We obtained tissues from Western Scrub-Jays represent- 
ing all known subspecies, genetic lineages, and geographic 
regions from existing frozen tissue specimens in museums 
or from new specimens (Table 2; see [30] for a complete 
list of the 689 specimens used in this study and associated 
data). We collected new specimens with United States fed- 
eral permit MB101618 and Nevada state permit S3505. 
Field collecting complied with Occidental College's Insti- 
tutional Animal Care and Use Committee protocol R12- 
1011-01. We then extracted DNA from blood or tissue 
using a DNeasy Tissue Kit (QIAGEN Inc., Valencia, CA) 
according to the manufacturer's protocols. 

Assigning coastal vs. interior mtDNA haplotypes 

Phylogenies based on mtDNA have been characterized in 
previous studies of Western Scrub-Jays using a limited 
numbers of individuals [23,24]. Our goal was to conduct 
deep mtDNA sampling across the range of the species and 
assign coastal and interior types to determine the extent 
of mixing within populations, and to determine whether 
rare long-distance dispersal events have occurred. 

Instead of direct sequencing, we screened individuals for 
coastal versus interior mtDNA type with a rapid, diagnos- 
tic assay. We used previously sequenced individuals [24] 
to identify a rare-cutting restriction enzyme that differed 
between coastal individuals, which had the cut site, and 
interior individuals, which did not. Then, for new samples 
of unknown haplotype, we amplified the cyt b gene using 
previously published primers and PCR conditions [24]. 
Once completed, we added 0.5 uL of BSR-DI (2000 U/mL, 
New England Biolabs, Ipswitch, MA) to each reaction and 
incubated for one hour at 65°C, followed by visualization 
on an agarose gel. We scored individuals having two bands 
as coastal and individuals having one band as interior, 
employing both positive and negative controls. We ran 
samples where mtDNA type conflicted with geographic 
(sampling) location a second time to confirm the result. 

Nuclear microsatellite loci 

We used microsatellites to gain insight into differenti- 
ation and gene flow of the nuclear genome and to score 
individuals as potential hybrids (given that mtDNA is 
haploid and provides no information as to the hybrid an- 
cestry of individuals). Previously characterized microsat- 
ellite loci have been isolated from the Florida Scrub-Jay 
[48]. From these loci, we chose a panel of 14 loci that 
amplified well and were variable among Western Scrub- 
Jays. To reduce processing time and cost while maintain- 
ing data quality, we pooled up to three loci and amplified 
them together (multiplexing) with a single dye [49], after 
first confirming that allele size classes for pooled loci did 
not overlap. 
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Table 2 Locality information for Western Scrub-Jays sampled for this study 



Map 


Location 


Coordinates 


n 


Subspecies 


I 


Vida, Oregon, USA 


44.1463, 


-122.5701 


6 


immamis 


2 


Alturas, California, USA 


41.55, -1 


20.6667 


23 


superciliosa 


3 


Lassen Co., California, USA 


40.7047, 


-120.7391 




superciliosa 


4 


Douglas City, California, USA 


40.5833, 


-123 


18 


caurina 


5 


Beckwourth, California, USA 


39.75, -1 


20.3667 




superciliosa 


6 


Poole Station Rd, Calaveras Co., California, USA 


38.1,-120.6333 


17 


superciliosa 


/ 


Martinez, Contra Costa Co., California, USA 


38.0128, 


-122.1325 




oocleptica 


8 


El Sobrante, Contra Costa Co., California, USA 


37.9755, 


-122.2859 


] 


oocleptica 


9 


Alhambra Valley Rd., Contra Costa Co., California, USA 


37.9663, 


-122.2059 


] 


oocleptica 


10 


Berkeley, California, USA 


37.8788, 


-122.2583 


] 


oocleptica 


1 1 


Carmel Valley Rd., Alameda Co., California, USA 


37.874, - 


121.832 




superciliosa 


12 


Moraga 1, Contra Costa Co., California, USA 


37.8598, 


-122.1343 




oocleptica 


12 


Moraga 2, Contra Costa Co., California, USA 


37.8653, 


-122.1386 




oocleptica 


12 


Moraga 3, Contra Costa Co., California, USA 


37.8602, 


-122.1235 




oocleptica 


12 


Orinda, Contra Costa Co., California, USA 


37.8654, 


-122.1519 


17 


oocleptica 


13 


Danville, Contra Costa, Co., California, USA 


37.8307, 


-121.9826 


1 


oocleptica 


14 


San Ramon, Contra Costa Co., California, USA 


37.75, -1 


22 


26 


oocleptica 


15 


Dublin, Alameda Co., California, USA 


37.7131, 


-121.9247 




oocleptica 


16 


Purisima Creek, San Mateo Co., California, USA 


37.4070, 


-122.4059 




oocleptica 


17 


North Fork, Madera Co, California, USA 


37.2429, 


-119.5421 




superciliosa 


18 


Independence, Inyo Co., California, USA 


36.7667, 


-118.2833 


6 


superciliosa 


19 


Bradley, Monterey Co., California, USA 


35.8690, 


-120.9267 


21 


californica 


20 


Bodfish, Kern Co., California, USA 


35.6167, 


-118.5 


20 


superciliosa 


21 


Big Bear City, San Bernadino Co, California, USA 


34.3167, 


-116.8333 


5 


obscura 


22 


Pioneertown, San Bernadino Co, California, USA 


34.1667, 


-116.5333 


2 


obscura 


23 


Santa Ana Mountains, Orange Co, California, USA 


33.7, -11 


7.6167 


16 


obscura 


24 


Love Valley, San Diego Co, California, USA 


33.5833, 


-116.7833 


12 


obscura 


25 


Pinyon Flat, Riverside Co, California, USA 


33.5833, 


-116.4667 


6 


obscura 


26 


Vallecito Mountains, San Bernadino Co, California, USA 


33.027,-1 


16.242 


10 


obscura 


27 


Mountain Springs Pass, Imperial Co, California, USA 


32.6667, 


-116.0833 


1 


obscura 


28 


Jacumba, San Diego Co, California, USA 


32.6333, 


-116.2167 


/ 


obscura 


29 


La Rumorosa, Baja California, Mexico 


32.55, -1 


16.05 


8 


obscura 


30 


Potrero, San Diego Co, California, USA 


32.65, -1 


16.6167 


3 


obscura 


31 


La Rosa de Castilla, Baja California, Mexico 


32.05, -1 


16.1333 


20 


obscura 


32 


San Lucas, Baja California Sur, Mexico 


27.5, -11 


2.3333 


1 


cactophila 


33 


Bahia Magdalena, Baja California Sur, Mexico 


24.7833, 


-112.1 


15 


cactophila 


34 


La Burrera, Baja California Sur, Mexico 


23.5, -11 


0.1167 


12 


hypoleuca 


35 


Virginia Mountains, Storey Co, Nevada, USA 


39.4977, 


-119.6344 


4 


superciliosa x nevadae 


36 


Lousetown, Storey Co, Nevada, USA 


39.3854, 


-119.6322 


1 


superciiiosa x nevadae 


37 


Brunswick Cyn, Pine Nut Mountains, Douglas Co, Nevada, USA 


39.1458, 


-119.6491 


4 


superciliosa x nevadae 


38 


Sunrise Pass Rd. 1, Pine Nut Mountains, Douglas Co, Nevada, USA 


39.066, - 


119.569 


2 


superciliosa x nevadae 


39 


Sunrise Pass Rd. 2, Pine Nut Mountains, Douglas Co, Nevada, USA 


39.0654, 


-119.5586 


2 


superciliosa x nevadae 


40 


Lebo Spring, Pine Nut Mountains, Douglas Co, Nevada, USA 


39.0557, 


-119.5703 


22 


superciliosa x nevadae 


41 


Sunrise Pass Rd. 3, Pine Nut Mountains, Douglas Co, Nevada, USA 


39.0523, 


-119.6028 


2 


superciliosa x nevadae 
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Table 2 Locality information for Western Scrub-Jays sampled for this study (Continued) 



42 


Hot Springs Mountain, Douglas Co., Nevada, USA 


39.0333, 


-119.6333 


6 


superciiiosa x nevadae 


43 


Pine Nut Creek, Douglas Co., Nevada, USA 


38.9414, 


-1 19.7069 


6 


superciliosa x nevadae 


44 


Gardnerville, Douglas Co., Nevada, USA 


38.8333, 


-1 19.6167 


5 


superciiiosa x nevadae 


45 


Woodford, Alpine Co, California, USA 


38.7833, 


-1 19.75 


/ 


superciliosa x nevadae 


46 


Nye Canyon, Pine Grove Hills, Lyon Co., Nevada, USA 


38.5703, 


-1 19.2083 


10 


nevadae 


4/ 


Pine Grove Hills, Lyon Co, Nevada, USA 


38.5133, 


-1 19.2040 




nevadae 


48 


Lucky Boy Pass Rd, Wassuk Range, Mineral Co, Nevada, USA 


38.4549, 


-118.6701 




nevadae 


49 


North Canyon, Wassuk Range, Mineral Co, Nevada, USA 


38.4625, 


-118.6475 




nevadae 


50 


Excelsior Mountains, Mineral Co, Nevada, USA 


38.2792, 


-118.4844 




nevadae 


51 


White Mountains, Inyo Co, California, USA 


37.2833, 


-118.1667 


16 


nevadae 


52 


Mt Charleston, Clark Co, Nevada, USA 


36.3667, 


-115.6333 


27 


nevadae 


53 


Toiyabe Mountains, Lander Co, Nevada, USA 


39.3333, 


-117.1333 


15 


nevadae 


54 


Baker, Millard Co, Utah, USA 


39.05, -1 


14.0833 




nevadae 


bb 


Hualapi Mountains, Mohave Co, Arizona, USA 


35.15, -1 


13.9 




nevadae 


56 


Drake, Yavapai Co, Arizona, USA 


35,-112.25 


26 


nevadae 


SI 


Stansbury Mountains, Tooele Co, Utah, USA 


40.35, -1 


12.5333 


15 


woodhouseii 


58 


Uinta Mountains, Duschesne Co, Utah, USA 


40.5816, 


-110.0105 




woodhouseii 


59 


Duschesne Co, Utah, USA 


39.9813, 


-110.2620 


1 


woodhouseii 


60 


Bonanza, Uitah Co, Utah, USA 


39.7811, 


-109.0260 


] 


woodhouseii 


61 


Whitewater, Mesa Co, Colorado, USA 


38.9199, 


-108.4834 




woodhouseii 


62 


Monument, Teller Co, Colorado, USA 


39.0938, 


-104.8296 


, 


woodhouseii 


63 


Gardner, Huerfano Co, Colorado, USA 


37.8833, 


-105.2 


13 


woodhouseii 


64 


Manzano, Valencia Co, New Mexico, USA 


34.6667, 


-106.4667 


26 


woodhouseii 


65 


Fort Davis, Jeff Davis Co, Texas, USA 


30.7, -104.1333 


22 


woodhouseii 


66 


Carta Valley, Edwards Co, Texas, USA 


29.8333, 


-100.6833 


28 


texana 


67 


Villa Ocampo, Durango, Mexico 


26.4667, 


-105.4833 


12 


grisea 


68 


Sombrerete, Zacatecas, Mexico 


23.7, -103.75 


9 


grisea 


69 


Rancho Santa Rita, Jalisco, Mexico 


21.45, -1 


01.9167 


19 


grisea 


70 


El Diamante Pass, Coahuila, Mexico 


25.3667, 


-100.8667 


21 


cyanotis 


71 


Bledos, San Luis Potosi, Mexico 


21.8667, 


-101.15 


20 


cyanotis 


72 


Taxco, Guerrero, Mexico 


18.5833, 


-99.6333 


3 


remota 


73 


Xocomanantlan, Guerrero, Mexico 


1 7.55, -99.65 


12 


remota 


74 


San Lorenzo de Abarrados, Oaxaca, Mexico 


17, -96.1667 


24 


sumichrasti 



Map numbers correspond to sampling sites in figures. 



After amplification, we determined fragment size using 
an ABI PRISM 3730 capillary sequencer and analyzed the 
resulting data using Geneious v. 6.0.4 (Biomatters). We 
called alleles objectively by creating bins based on known 
microsatellite repeat motifs [48], and assessed artifacts 
after genotyping (e.g., stuttering, null alleles, big-allele 
drop-out) using Microchecker [50]. To ensure that micro- 
satellite loci were evolving under neutral processes, we 
assessed each locus for deviation from Hardy- Weinberg 
equilibrium (HWE) using GenePop on the Web (http:// 
genepop.curtin.edu.au). Strong population structure can 
lead to false positives for deviations from HWE when 



populations are pooled (Wahlund Effect), so we assessed 
HWE in each population at each locus individually, apply- 
ing Bonferroni correction to control for false positives 
resulting from the effect of conducting multiple simultan- 
eous tests. We also performed a test of linkage disequilib- 
rium to confirm independent sorting of loci. 

Analysis of microsatellite variation for genetic clustering 
across the geographic range 

We assessed genetic structure of microsatellite data using 
Structure v. 2.3.4 [51], which uses Bayesian analysis to 
infer the number of populations (K) from a group of 
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individuals based only on their genetic variation, without 
prior information on where the individuals originated. We 
did not use the method of Evanno et al. [52] to detect the 
'true K because this method is known to underestimate 
genetic structure in all but cases of very strong genetic dif- 
ferentiation [53]. Initial tests of the Evanno et al. [52] 
method suggested K = 2 as the best model for our data; 
however, analysis at higher K revealed strong structure 
that was highly justified based on biology and geography. 
Thus, because one of our goals was discovery of the basic 
units of spatially structured genetic diversity, we instead 
used Structure to analyze successively smaller clusters of 
individuals at K = 2 until Structure could no longer find 
any clustering within groups, as demonstrated previously 
[54]. For example, when two clear geographic population 
clusters with strong biological justification emerged in the 
complete data set, we ran each of these clusters in its own 
analysis until K = 2 revealed no remaining geographical 
clusters. After the first run, we excluded hybrid popula- 
tions (identified as populations with both coastal and in- 
terior mtDNA types) from the analysis. We ran Structure 
with an admixture model, correlated allele frequencies, 
and a burn-in period set to 100,000 generations, followed 
by 500,000 generations, which was sufficient for each run 
to reach stationarity. 

Analysis of microsatellite variation across the hybrid zone 

To assess population assignment and potential hybrid an- 
cestry of individuals across the hybrid zone, we delimited 
13 populations along a transect running northwest to 
southeast from the Sierra Nevada in northeastern California 
into west-central Nevada (Figure 2). This line conforms to 
prior descriptions of the contact zone [25,55,56] and our 
own field observations. To obtain estimates of hybrid an- 
cestry (i.e., Q-scores) within the hybrid zone, individuals 
from these 13 populations were analyzed in Structure sep- 
arately from all other individuals in the study using K = 2; 
their resulting Q-scores were recorded for later use in cline 
analysis. The Q-score values provide an overall estimate of 
nuclear variation that is expected to suffer less from ran- 
dom processes affecting individual loci. Moreover, cline 
analysis requires binary data or frequency values varying 
from 0 to 1. Due to large number of alleles at a given locus, 
individual microsatellite loci do not usually conform to this 
standard. 

Testing for phylogenetic structure of microsatellite 
groups using mtDNA 

Previous studies that generated mtDNA phylogenies of 
Western Scrub-Jays did not sample deeply within popula- 
tions [23,24]. Thus, to determine if nuclear DNA groups 
discovered in this study were discernible using mtDNA, 
we sequenced multiple individuals from all major groups 
suggested by our microsatellite results for the mtDNA cyt 



b gene, using previous protocols [24]. We combined these 
new data with existing data for Aphelocoma and close out- 
groups {Calocitta formosa, Gymnorhinus cyanocephalus, 
and Cyanocitta stelleri) from GenBank, and we generated 
a phylogeny from these sequences using BEAST v. 1.7.5 
[57]. We ran the analysis using the best-fit substitution 
model for the cyt b gene determined through Akaike In- 
formation Criterion (AIC), which was HKY + I + G with 
six gamma categories and base frequencies estimated as 
determined by MrModeltest, v2 [58]. We used a strict 
clock with a random starting tree. We ran the chain to 
1.0 x 10 6 iterations, sampling every 1,000. The analysis 
converged quickly as determined by observing the plot 
of likelihood scores and ESS scores > 200 with Tracer. 
We discarded the first 1,000 of 10,000 trees as burn-in and 
produced a consensus tree from the remaining posterior 
trees using Geneious [59]. The maximum clade credibility 
produced with TreeAnnotator was almost identical, 
with key nodes mentioned in the Results having the 
same PP values. 

Fitting geographic dines to genetic data across the 
hybrid zone 

We used HZAR v2.5 [60], a statistical package imple- 
mented in R, to determine the width of geographic 
clines estimated from mtDNA and nuclear markers, and 
the extent of mismatch, if any. We measured the dis- 
tance between points from the coastal northwest end of 
the hybrid zone (site 2 in Figure 1), through the contact 
zone, continuing to pure interior populations to the 
southeast (site 52 in Figure 1). These were then com- 
pressed to a single line using HZAR. This transect likely 
captures the major axis of gene flow across the hybrid 
zone [25], with minor caveats covered in the Discussion. 
We estimated clines for mtDNA haplotype frequency 
and average Q-score per site as determined in our Struc- 
ture analyses. We used these cline analyses to estimate 
changes of the molecular characters in local mean fre- 
quencies. We modeled the cline shape using three equa- 
tions [61,62] describing a sigmoid shape at the center of 
the transition with two exponential decay curves on ei- 
ther side of the transition. We estimated several parame- 
ters, including width (w), center (c), delta (d, distance 
between the center and the tail), and tau (t, slope of the 
tail). We also incorporated the possibility that Pmax and 
Pmin (the "top" and "bottom" of the cline) were either 
fixed or free to vary. We fit three sets of five cline 
models using the Metropolis-Hasting algorithm in R. 
Model set 1 has no scaling {Pmin = 0, Pmax =1), model 
set 2 has fixed scaling {Pmin = observed minimum, 
Pmax = observed maximum), and model set 3 allows 
Pmin and Pmax to vary. Within each model set, scaling 
and tails are fixed or free to vary. We compared these 
models to a null model of no clinal transition. 
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For the two genetic datasets, we set up a covariance 
matrix by running each model for 1.0 x 10 6 generations. 
We ran three independent chains for 9.0 x 10 6 generations 
and assessed for convergence and stability by visualizing 
the MCMC traces. We performed model selection using 
corrected AIC (AICc). We estimated two log-likelihood 
confidence intervals around each parameter estimate. If 
the parameter estimate from one dataset did not overlap 
with the confidence interval of the parameter estimate 
from the other dataset, we considered them to be signifi- 
cantly different. We also used estimates of cline width to 
determine the time since contact under a neutral diffusion 
model solving for t in the equation w = 2.51cn/i [63], where 
t is the time since contact and a is the root mean square 
dispersal distance from estimates in [31]. Because females 
and males have different dispersal distances, we used the 
female dispersal distance (1.5 - 5 km) to estimate time 
since contact using mtDNA, and we used the average dis- 
persal estimate of males and females (0.5 - 3.0) to esti- 
mate time since contact for the nuclear dataset. 

Correlating phenotype with hybrid ancestry 

We collected phenotypic data on a subset of 124 coastal 
and interior Western Scrub-Jay specimens within and 
outside the hybrid zone. Unfortunately we could not 
perform a cline analysis along the transect because we 
did not have sufficient sample sizes from pure popula- 
tions at both ends of the transect. However, to deter- 
mine the extent to which phenotype is correlated with 
nuclear DNA variation, we conducted regression analysis 
of phenotypic traits (plumage and morphology) with Q- 
scores at several spatial scales. The set of 124 specimens 
includes the full ranges of coastal and interior pheno- 
types, many hybrids, and pure individuals of each lineage 
with no evidence for hybrid genetic ancestry (Q-score at 
or close to 1 or 0). Because environment can influence 
phenotype when considering allopatric lineages, a more 
robust test of the genotype-phenotype link occurs within 
the hybrid zone where individuals share the same envir- 
onment and there is more continuous variation in Q- 
scores. Here, we looked specifically at the 68 specimens 
that are part of the hybrid transect, as well as a more re- 
strictive sample of 50 individuals within the core hybrid 
zone defined as sites 35-45 in Figure 1. 

To ensure consistent measurements, one author (TCW) 
measured wing, tail, tarsus, bill length, bill depth, and bill 
width to the nearest 0.1 mm with digital calipers on all 
124 birds after first verifying high repeatability scores for 
all traits. On a subset of 66 of these birds, another author 
(JEM) assessed three qualitative plumage traits known to 
vary between coastal and interior lineages: the amount of 
blue collar (reduced in interior birds); size of the eyestripe 
(reduced in interior birds); and color of the undertail co- 
verts (blue-tinged in interior birds as opposed to white in 



coastal birds). We determined variation in plumage traits 
by lining up all 66 specimens at the same time and arran- 
ging them in order of trait appearance (e.g., very white 
undertail coverts to very blue undertail coverts). We then 
divided the arranged specimens into six even categories 
and assigned a value of 1 to 6. We repeated these steps for 
each trait. We analyzed morphological data with principal 
components analysis on the correlation matrix in Stata 10 
to identify axes defining the most variation in the data. 
We assessed relationships between Q-scores and univari- 
ate and multivariate phenotypic traits with linear regres- 
sion at each of the three spatial extents described above 
(full range, hybrid transect, core hybrid zone) using Stata. 

Availability of supporting data 

New DNA sequences have been deposited in Genbank 
under accession numbers KJ835799-KJ835861. Museum 
catalog numbers, localities, raw microsatellite data, pheno- 
typic data, and Q-scores for each individual are available 
through Dryad http://dx.doi.org/10.5061/dryad.57f48. 
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Additional file 1: Loadings for principal components analysis. 
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